#!/usr/bin/python
# GM_N_bvals.py
'''Condenses N to relevant points by linear interpolation or extrapolation.

In the case of an altered MvList (if interpolateMV is True), also sums over Mv.'''

import math
import numpy

import GM_InterpRoutines as IR

VariableNames = ['ba', 'vanillalimitset', 'limitset', 'N_bM', 'bvals']

def N_bvals(simpars, phypars, N):
    # Bear in mind that 'N' is log10 the stars per deg**2
    # Input
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    N = 10**N
    # Relevant distances
    ba = ba*(180./math.pi) # radians -> degrees
    # Sum over magnitudes if necessary
    N = IR.DownSampleLinear(vanillalimitset, limitset, N)

    # Interpolate or extrapolate
    compN = numpy.zeros(N_bM.shape, dtype=float)
    for i in range(0, len(bvals)):
        # Find where ba ~ bvals[i]
        b = bvals[i]
        mask1 = numpy.less(ba, b)
        n0 = numpy.array(mask1, dtype=int).sum()
        if n0 != 0:
            n0 = n0 - 1
        if n0 == (len(ba) - 1):
            n0 = n0 - 1
        n1 = n0 + 1
        N0 = N[:,n0] # 1d
        N1 = N[:,n1] # 1d
        b0 = ba[n0] # 0d
        b1 = ba[n1] # 0d
        (slope, offset) = IR.linear(b0,N0,b1,N1)
        thisN = slope*b + offset

#        db0 = b0 - b # 0d
#        db1 = b1 - b # 0d
#        N0db1 = N0*db1 # 1d
#        N1db0 = N1*db0 # 1d
#        num = N0db1 - N1db0 # 1d
#        denom = db1 - db0 # 0d
#        thisN = num/denom
        compN[:,i] = thisN
    compN = numpy.log10(compN)
    return compN
